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ABSTRACT 

We present a Monte Carlo radiative transfer technique for calculating synthetic 
spectropolarimetry for multi-dimensional supernova explosion models. The approach 
utilises “virtual-packets” that are generated during the propagation of the Monte 
Carlo quanta and used to compute synthetic observables for specific observer orien¬ 
tations. Compared to extracting synthetic observables by direct binning of emergent 
Monte Carlo quanta, this virtual-packet approach leads to a substantial reduction 
in the Monte Carlo noise. This is vital for calculating synthetic spectropolarimetry 
(since the degree of polarisation is typically very small) but also useful for calcula¬ 
tions of light curves and spectra. We first validate our approach via application of an 
idealised test code to simple geometries. We then describe its implementation in the 
Monte Carlo radiative transfer code ARTIS and present test calculations for simple 
models for Type la supernovae. Specihcahy, we use the well-known one-dimensional 
W7 model to verify that our scheme can accurately recover zero polarisation from a 
spherical model, and to demonstrate the reduction in Monte Carlo noise compared 
to a simple packet-binning approach. To investigate the impact of aspherical ejecta 
on the polarisation spectra, we then use ARTIS to calculate synthetic observables for 
prolate and oblate ellipsoidal models with Type la supernova compositions. 

Key words: polarisation - radiative transfer - methods: numerical - super novae: 
general 


1 INTRODUCTION 

Type la supernovae (SNe la) are generally believed to be 
thermonuclear explosions of carbon-oxygen white dwarfs 
(see e.g. Ropke et al. 2011 and Hillebrandt et al. 2013 for 
reviews). However, answers to the questions of how and why 
the explosion is triggered remain unclear. Most of the estab¬ 
lished theoretical models involve close binary systems, but 
we still do not know if the companion star is a second white 
dwarf (double degenerate system, Webbink 1984; Iben Sz 
Tutukov 1984) or a non-degenerate star (single degenerate 
system, Whelan & Iben 1973), and whether the explosion 
is triggered when an accreting white dwarf approaches the 
Chandrasekhar limit or via some other process. 

For Chandrasekhar-mass models, neither a pure defla¬ 
gration nor a pure detonation model is able to fully account 
for the properties observed in SNe la: the former leads to 
strong turbulence and buoyancy resulting in fingers of nickel 
and carbon-oxygen at all ejecta velocities (Gamezo et al. 
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2003; Ropke et al. 2006; Jordan IV et al. 2012; Ma et al. 
2013; Fink et al. 2014) while the latter fails to produce 
intermediate-mass elements (Arnett 1969). However an in¬ 
terplay of these two models, so-called delayed-detonation 
models, remains promising for providing a good match to 
data. Possibilities include spontaneous deflagration to det¬ 
onation transition models (Khokhlov 1991; Hoflich et al. 
1995, 1996; Hoflich V Khokhlov 1996; Kasen et al. 2009; 
Blondin et al. 2012), and the gravitationally conhned det¬ 
onation model (Plewa et al. 2004; Jordan IV et al. 2008, 
2012). For non-Chandrasekhar-mass models, in which the 
accreting white dwarf may explode without approaching the 
Chandrasekhar limit, viable mechanisms are the detonation 
of helium layers on the surface of the accreting white dwarf 
(the double detonation model; see Nomoto 1980; Woosley 
et al. 1980; Livne 1990; Woosley V Weaver 1994; Fink et al. 
2007, 2010; Shen V Bildsten 2009; Woosley V Kasen 2011; 
Moll V Woosley 2013) and the violent mergers of two white 
dwarfs (Pakmor et al. 2010, 2012, 2013; Moll et al. 2014; 
Raskin et al. 2014). Alternative models include the explosion 
of white dwarf merger remnants (Benz et al. 1990; Van Kerk- 
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wijk et al. 2010; Shen et al. 2012; Kashyap et al. 2015), the 
head-on collisions of white dwarfs (Rosswog et al. 2009), pos¬ 
sibly induced in triple systems (Kushnir et al. 2013), or the 
merger of a white dwarf with the hot core of an asymptotic 
giant branch star (Soker et al. 2014). 

There are a variety of ways to attempt to determine 
which of the proposed progenitor/explosion channels really 
occur (see e.g. Maoz et al. 2014 for a review). One approach 
is to perform explosion simulations with associated radia¬ 
tive transfer calculations and compare their predictions to 
data. Such work is well established (e.g. Hoflich et al. 1993; 
Hauschildt & Baron 1999; Kasen et al. 2006; Kromer & Sim 
2009; Blondin et al. 2012; Dessart et al. 2014; Wollaeger 
& Van Rossum 2014) and it is now possible to compute 
synthetic observables from multi-dimensional models for a 
variety of explosion scenarios. Although state-of-the-art ex¬ 
plosion simulations allow us to capture considerable com¬ 
plexity in SN la models (e.g. turbulence), degeneracies be¬ 
tween models remain and make unambiguous interpretation 
difficult, even for the best observed nearby examples (Ropke 
et al. 2012). 

One potentially powerful discriminant is the geome¬ 
try, which can be quite different between models and de¬ 
pend on the nature of the progenitor and explosion mecha¬ 
nism. From the observational side, both nebular phase spec¬ 
troscopy (Gerardy et al. 2007; Maeda et al. 2010) and spec- 
tropolarimetric observations (see Wang V Wheeler 2008 for 
a review) provide evidence that SNe la are not perfectly 
spherically symmetric. Continuum polarisation is typically 
quite low (0.2 — 0.3 per cent) in normal SNe la prior to op¬ 
tical maximum, pointing toward very small departures from 
global spherical symmetry; significant polarisation is, how¬ 
ever, found across the line profiles of spectral features as¬ 
sociated with intermediate-mass elements (calcium, silicon, 
sulphur and magnesium, but not oxygen; see e.g. Wang V 
Wheeler 2008 and reference therein), suggesting that asym¬ 
metries in the element distribution are present. Other sub¬ 
classes of SNe la seem to display some peculiarities: sub- 
luminous SNe la show higher continuum polarisation levels 
(0.3 — 0.8 per cent, Howell et al. 2001; Patat et al. 2012), 
whereas high-velocity SNe la show stronger line polarisation 
(~ 2 per cent, Leonard et al. 2005; Wang et al. 2006). 

Asymmetric ejecta for SNe la are also predicted by 
many multi-dimensional explosion models, although with 
different degrees and types of asymmetry. For instance, the 
delayed-detonation model predicts that the ejecta can be 
quasi-spherical on large scales but with complex substruc¬ 
tures (in both density and composition, Seitenzahl et al. 
2013; Sim et al. 2013), whereas simulations of violent white 
dwarf mergers predict departures from spherical symmetry 
on large angular scales (Pakmor et al. 2010; Moll et al. 2014; 
Raskin et al. 2014). These differences in ejecta geometry be¬ 
tween models have led to suggested connections with ob¬ 
served objects. For example, Patat et al. (2012) proposed 
that the sub-luminous SN 2005ke might be explained by an 
explosion of a rotating white dwarf or a double-degenerate 
merger. It has also been suggested that several models could 
be ruled out because they yield explosions that are too as- 
pherical and therefore inconsistent with continuum polarisa¬ 
tion measurements for normal SNe la. For instance, Maund 
et al. (2013) have claimed that the low continuum polarisa¬ 
tion and significant line polarisation observed in SN 2012fr is 


consistent with delayed-detonation models but inconsistent 
with deflagration models or violent white dwarf mergers. 

However, interpreting polarisation data and quantify¬ 
ing arguments about its implications for models is difficult 
because estimating the degree of polarisation expected for 
a complex ejecta morphology (as provided e.g. by multi¬ 
dimensional explosion simulations) is not trivial: polarisa¬ 
tion depends on the opacity distributions in a complex way 
(Hoflich 1991; Dessart V Hillier 2011). To date, spectropo- 
larimetric data of SNe la have often been interpreted by 
comparing the observed polarisation levels with predictions 
from toy models with idealised configurations for the ejecta, 
e.g. ellipsoidal and clumped shell models or spherical shells 
with a hole or toroid (Howell et al. 2001; Kasen et al. 2003, 
2004; Patat et al. 2012). These idealised geometries are 
well-suited for building intuition and establishing the frame¬ 
work for interpreting polarisation data. However, quantita¬ 
tive comparisons between predictions of multi-dimensional 
explosion models and data require that polarisation calcula¬ 
tions are made for the complete density/composition distri¬ 
butions predicted by hydrodynamic simulations. Such calcu¬ 
lations make it possible to quantitatively compare the pre¬ 
dictions of models to data (and each other) and assess the 
extent to which their geometries can really be distinguished 
via polarisation. 

Here we present a polarisation scheme recently im¬ 
plemented in the three-dimensional, time-dependent Monte 
Carlo radiative transfer code ARTIS (Applied Radiative 
Transfer In Supernovae, Sim 2007; Kromer & Sim 2009). The 
scheme involves two parts, first an implementation of Stokes 
parameters for the Monte Carlo quanta and secondly the de¬ 
velopment of techniques to reduce the Monte Carlo noise in 
the emergent synthetic observables. This is particularly im¬ 
portant when we aim to extract very weak polarisation sig¬ 
nals (low percentage levels are observed in SNe la), but also 
useful for total flux spectra and light curves. Although in 
this work we focus on the development, implementation and 
validation of the method using one-dimensional and two- 
dimensional models, our particular technique is well suited 
to exploit the multi-dimensional capability of ARTIS and 
therefore to be applied to three-dimensional explosion mod¬ 
els. Details of the methodology used are given in Section 
2 and the polarisation scheme is validated via an idealised 
test code in Section 3. We then present first results from 
the implementation in ARTIS, including testing with the 
one-dimensional W7 model (Nomoto, Thielemann & Yokoi 
1984; Iwamoto et ah, 1999) and two-dimensional ellipsoidal 
toy models (Section 4). We summarise and draw conclusions 
in Section 5. 


2 METHOD 

The polarisation of a beam of radiation is characterised by 
the four-dimensional Stokes vector S — The 

first component, /, gives the total intensity, Q and U mea¬ 
sure the degree of linear polarisation and V of circular polar¬ 
isation. Since circular polarisation has never been observed 
in SNe la, and the radiative transfer calculations for circular 
and linear polarisation can be decoupled in scattering at¬ 
mosphere in the absence of magnetic fields (Chandrasekhar 
1960), here we neglect the V component. The Stokes vector 
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Figure 1. The meridian plane coordinate system adopted in the 
Monte Carlo code. The Stokes vector is defined in the plane or¬ 
thogonal to the direction of propagation, n. Q is defined as the 
intensity difference between two perpendicular reference axes, I 
and r, whereas U is the equivalent difference with I and r counter¬ 
clockwise rotated by 45 degrees (viewed antiparallel to n). 


is defined in the plane orthogonal to the direction n in which 
the radiation propagates. To define the Stokes parameters 
for linearly polarised radiation, we introduce two reference 
axes I and r so that I lies in the meridian plane (plane de¬ 
fined by n and the polar axis z) and n — r X I (see Fig. 
1). With this convention, Q is defined as the difference be¬ 
tween intensity I\ with electric field oscillating along I and 
intensity Ir with electric field oscillating along r; U is the 
equivalent difference in intensities with the reference axes I 
and r counter-clockwise rotated by 45 degrees (as viewed 
looking antiparallel to n) to give a and b. The resulting 
Stokes vector S can be expressed as 





//l + /r 
/l - /r 

y/a - h 






( 1 ) 


or in terms of a dimensionless Stokes vector, s: 



The polarisation fraction, p, and the position angle, of a 
beam are related to the Stokes parameters by 


2.1 Propagation 

In this section, we adopt the terminology introduced by 
Lucy (2002, 2005) and discuss the general scheme used to 
include polarisation in our radiative transfer code. The cal¬ 
culations we present use Monte Carlo methods: the radiative 
transfer problem is solved by simulating the propagation of 
Monte Carlo quanta (packets of identical photons) through 
an expanding medium. The propagation of an r—packet 
(monochromatic packet of ultraviolet-optical-infrared radia¬ 
tion) is followed through the ejecta in the rest frame (rf) and 
stopped by interactions with matter (which are treated in 
the comoving frame, cmf). For ultraviolet-optical-infrared 
radiation, the code currently accounts for both line opac¬ 
ity (treated in the Sobolev approximation [Sobolev I960]) 
and continuum opacity due to electron scattering, bound- 
free and free-free absorption. To choose when continuum 
and interaction events occur we use the method outlined 
by Mazzali & Lucy (1993). Once a random optical depth Tr 
is drawn, we determine the trajectory point at which the 
r—packet interacts with the next line (using the Sobolev 
approximation): if the continuum opacity accumulated up 
to that point is greater than Tr, a continuum absorption 
is selected; if instead the sum of continuum and line opac¬ 
ity is greater than Tr, a line event occurs; otherwise, the 
process is repeated for the next line with which the packet 
comes into resonance. In an electron scattering event, the 
r—packet keeps the same cmf frequency and is assigned a 
new direction of propagation. For all other interactions (line 
absorption, free-free absorption and bound-free absorption) 
either a /c—packet (packet of thermal kinetic energy) or an 
z—packet (packet of excitation/ionisation energy) is acti¬ 
vated and then processed according to the scheme proposed 
by Lucy (2002) as described by Kromer & Sim (2009). 

Our polarisation scheme adopts a method similar to 
that proposed by Lucy (2005) and already implemented by 
Kasen et al. (2006). Polarisation is introduced by assigning a 
Stokes vector to each r—packet. When an interaction occurs, 
the Stokes parameters are transformed through the following 
sequence of steps: first, we transform the incoming Stokes 
vector Si in the rf to s( in the cmf ^. The plane in which the 
electric field oscillates changes as a result of the aberration 
of the direction n to 



,3) 

= 

where x is the angle between the electric field orientation 
and the reference axis 1. Spherically symmetric geometries 
are characterised by null polarisation since every contribu¬ 
tion is canceled by an orthogonal contribution one quadrant 
away, whereas aspherical geometries may lead to a polari¬ 
sation signal due to non perfect cancellation of the Stokes 
vectors (see also Kasen et al. 2003 and discussion in Section 

4). 


where c is the speed of light, 7 the Lorentz factor and v the 
local velocity of the ejecta (Castor 1972). To derive how the 
Stokes parameters change under the Lorentz transformation, 
we introduce a unit vector e that describes the orientation 
of the net electric field in the rf, 

e = (cos x) I - (sin x) r , ( 6 ) 

with the angle x between e and I computed from the incom¬ 
ing Stokes vector following equation (4). From this, it is easy 
to obtain a cmf Stokes vector representation that is relative 


^ In the following, cmf quantities are denoted with a prime, 
whereas unprimed quantities refer to the rf. 
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to axes I' and r' (defined hy n'), making use of the fact that 
the polarisation, p, is invariant (Cocke & Holm 1972): 

q[ = Pi cos 2x u{ = Pi sin 2%' • (7) 


Following transformation to the cmf, the Stokes param¬ 
eters can be updated in accordance with the physical inter¬ 
action that occurs. For hound-hound^ hound-free and free- 
free absorptions, we first activate either an z—packet or a 
/c—packet, as outlined in Kromer & Sim (2009) using the 
machinery described by Lucy (2002), and then convert this 
to a new r—packet. For all these processes, the r—packet 
is assumed to retain no information on polarisation and is 
reemitted in a random direction with zero polarisation^: 

q'f^O zzf = 0 . (8) 


If instead the r—packet undergoes eleetron seattering, we 
follow the scheme introduced by Chandrasekhar (1960) and 
discussed in terms of a Monte Carlo implementation by Code 
& Whitney (1995) and Whitney (2011). A scattering angle 
© is properly sampled (see below), a new direction nj is 
computed and the Stokes vector is transformed via the scat¬ 
tering matrix 

^ /cos^© + l cos^© — 1 0 \ 

A(©) = - j cos^© — 1 cos^© + 1 0 j . (9) 

y 0 0 2 cos©/ 


After applying the scattering matrix, the dimensionless 
Stokes vector is normalised so that its first component is 
equal to 1. As shown in Fig. 2, since A(©) is defined in the 
scattering plane, the Stokes vector is first rotated into this 
plane and then rotated back to the meridian frame after the 
scattering matrix A(©) has been applied, i.e. 


Sf = R{7r — 


where 




i2)A{e)R(- 


(10) 

0 

0 \ 


cos 20 

sin 20 

(11) 

— sin 20 

cos 20 / 



is the matrix to rotate the Stokes vector by an angle <f> clock¬ 
wise when facing the source. The rotation angles ii and 12 
are computed with the convention that q' — 1 represents an 
electric field oscillating in the scattering plane. The scatter¬ 
ing angle © is chosen by sampling the probability distribu¬ 
tion 


P(©,zi) = J |^(cos^©+l) + (cos^© —1)(^( cos 2zi —sin 2zi)j 

( 12 ) 

from equation (10) and using a reiection technique (Code & 
Whitney 1995). 

Finally the Lorentz transformation procedure (see 
above) with v ^ —v is used to transform the Stokes vector 
Si to S'f. 


^ Hamilton (1947) have shown that the angular distribution in 
resonant line scattering can be expressed as the sum of a dipole 
and an isotropic term (with relative contribution depending on 
the quantum numbers of upper and lower state of the line). Al¬ 
though polarisation may arise in cases where the dipole is the 
dominant term, its magnitude is typically lower than that from 
electron scattering and thus can be regarded as a second-order 
effect (Jeffery 1989, 1991). 


2; 



Figure 2. The geometry adopted for electron scattering in the 
coordinate system introduced by Chandrasekhar (1960) with the 
additional reference axes l[ and 1^ defined in the text. A packet 
moving along is scattered through an angle 0 to a new direc¬ 
tion nf The reference axis l[ (and accordingly the Stokes vector) 
are rotated through an angle ii into the scattering plane (in grey) 
and through an angle Z 2 after scattering. 


2.2 Spectrum extraction techniques 

Monte Carlo methods have been exploited in many radiative 
transfer calculations. These methods have proven particu¬ 
larly useful for multi-dimensional problems since the Monte 
Carlo algorithm can be easily implemented for arbitrary ge¬ 
ometry and scales very well for use on massively parallel 
compute systems. The main drawback of Monte Carlo meth¬ 
ods is their stochastic nature, which leads to solutions that 
are affected by Monte Carlo noise. Therefore, it is important 
to attempt to optimise Monte Carlo methods, particularly 
when we want to extract very weak signals (e.g. polarisa¬ 
tion) from the simulations. In this study we compare three 
different methods for extracting synthetic observables from 
our simulations, with the aim of selecting the most suitable 
to synthesise spectra with low Monte Carlo noise. In Section 
2.2.1 we present a simple direct counting approach, while in 
Section 2.2.2 and 2.2.3 we explore alternative techniques in¬ 
spired by Lucy (1999, 2005) and recently implemented by 
Kerzendorf & Sim (2014). 


2.2.1 Direet eounting teehnique 

In the Direct Counting Technique (DCT), Monte Carlo 
quanta are followed along their trajectories and their Stokes 
parameters updated at each interaction, following the pro¬ 
cedure outlined in Section 2.1. Packets that reach the outer 
boundary are collected in bins according to their final direc¬ 
tion n and frequency u and the resulting spectra are com¬ 
puted as 

^ Ai Ai/ 47rr2 

where Sf is the dimensionless Stokes vector of the escaping 
r—packet, e its rf energy, r the distance of the observer from 
the system and the sum is performed over all the r—packets 
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Figure 3. Sketches of the principle behind the EBT (left panel) 
and the TBT (right panel). EBT: for every interaction of the 
r—packet, a r’—packet is created and sent directly to a specific ob¬ 
server. Its Stokes parameters are treated in accordance with the 
specific interaction and are added as contributions to the emer¬ 
gent spectrum. TBT: contributions in the TBT come not only 
from trajectories that are terminated by physical interactions of 
the r—packet (as in the EBT, black points) but also from those 
terminated by numerical events (r—packets crossing boundaries, 
white points). The arc segment represents the outer boundary of 
the computational domain. 


that escaped in the selected angular bin, time interval 
[t —At/2, t+At/2] and frequency range [z/ —Az//2, z/+Az//2]. 

The direct counting approach provides a simple way of 
computing polarisation spectra for different viewing angles, 
i.e. for different observers. However, the need to average con¬ 
tributions from packets that escape in the same angular bin 
but with different angles inevitably leads to an approximate 
result: if the number of angular bins is too small, summing 
contributions from different angles will produce a poor esti¬ 
mate of the observables seen by a single observer; if instead 
too many angular bins are used, the number of packets es¬ 
caping per bin becomes small, leading to high Monte Carlo 
noise. 


2.2.2 Event based teehnique 

In this approach (Event-Based Technique, EBT) we still fol¬ 
low the propagation of Monte Carlo packets exactly as be¬ 
fore. However, whenever an r—packet interaction occurs the 
propagation is suspended, and a “virtual” packet (r^—packet, 
Kerzendorf & Sim 2014) is created and handled as described 
below. Once the i’—packet calculation is completed, the 
propagation of the original r—packet is resumed and this 
process repeated for every following interaction. This ap¬ 
proach is similar to, and inspired by, that used by Knigge 
et al. (1995) and Long & Knigge (2002). 

When a packet is created, it is always launched in the 
direction of a selected observer, riobs, and has cmf frequency 
and energy set equal to those of the r—packet at the point 
of creation (see Eig. 3). Since the r'—packet is forced to go 
to the observer, the scattering angle is now not determined 
by sampling either an isotropic distribution (z—packet or 
/c—packet deactivation) or equation (12) (electron scatter¬ 
ing), but rather calculated as 

cos 0 = n'n • n'bs . (14) 


The Stokes parameters, initially set equal to those of the 
incoming r—packet, are transformed in accordance with the 
physical process selected for the r—packet. If the r—packet 
is scattered by an electron, the cmf Stokes vector is trans¬ 
formed according to equation (10), with the scattering angle 
© determined by equation (14). If instead a new r—packet 
is created from an z—packet or /c—packet deactivation, we 
create an unpolarised z;—packet. 

Once created and assigned a rf frequency, a rf energy 
(e), a direction of propagation and a Stokes vector (sf), 
the z;—packet is propagated through the ejecta towards the 
observer and it is interpreted as a contribution to a bin 
in the emergent spectrum. Specifically, the flux (/) and 
polarisation (Q and U) spectra in a given frequency bin 
[jy — Az//2, jy + Aiyl2] and time interval [t — Atl2, t + At/2] 
can be expressed as 



w_^_ 

At Aiy 



(15) 


where the sum is performed over all the z;—packets escaping 
in the selected frequency and time bins, r is the distance 
of the observer from the system and the two terms inside 
parentheses are weighting factors accounting for the proba¬ 
bility that the z;—packet could reach the observer. Because 
the z;—packet is forced to go to the observer, we first account 
for the probability per unit solid angle associated with the 
chosen direction, that is: 


I ^ (fc-/i-)packet dead. 

l^P(0 ,ii) eledron scattering 

(16) 

The exponential factor in equation (15) accounts for the 
probability that the z;—packet could reach the observer 
without further interactions, with the optical depth to the 
boundary, Tesc, computed as 


Tesc 


Tcont 




"^sob 


(17) 


Here the sum is performed over all the line opacities (xsob) 
encountered by the z;—packet on its trajectory to the bound¬ 
ary and the continuum opacity, Tcont, is computed as an in¬ 
tegral of the continuum attenuation coefficient, /ccont, over 
trajectory length: 

Tcont — / ^cont ds . (18) 


The advantage of the z;—packet technique is that it al¬ 
lows us to compute spectra and light curves for any specific 
viewing angle, avoiding the need to average contributions 
from different angles in the same bin. Unlike in the direct 
counting approach, where an r—packet makes a single con¬ 
tribution to the emergent spectrum, the final spectra in the 
EBT also contain (appropriately weighted) information de¬ 
rived from every interaction that the r—packet undergoes 
(see Eig. 3). Eor these reasons, we expect this technique 
to produce spectra and light curves with lower Monte Carlo 
noise. However, the need of creating and handling z;—packets 
introduces a computational overhead that could make the 
EBT less efficient than the DCT. Such possibilities are ex¬ 
plored quantitatively in Section 3.2. 


© 0000 RAS, MNRAS 000 , 000-000 









































6 M. Bulla et al. 


2.2.3 Trajectory based technique 


In the EBT described above, we used all the r—packets to 
provide us with an ensemble of physical events, and then 
for each event in the ensemble we computed the probabil¬ 
ity of it giving rise to a photon escaping to the observer. In 
the third technique (Trajectory-Based Technique, TBT) we 
instead obtain from the r—packets an ensemble of photon 
trajectories, which can be taken as a discrete representation 
of the radiation held in the simulation. We can then estimate 
observables by summing over this ensemble of trajectories, 
computing for each the probability that interactions of radi¬ 
ation on that path could have given rise to photons escaping 
to the observer. 

This summation is achieved by generating a r’—packet 
for each trajectory path A/ (including those terminated 
by numerical events, e.g. packets crossing grid cell bound¬ 
aries) and sending it towards the observer at riobs 
(see Fig. 3). The i’—packet contribution to the emergent 
spectrum is hrst weighted by the probability per unit 
solid angle {dP/dTt\TBT) that photons on the path A/ 
could have undergone an interaction that gave rise to re¬ 
emission/scattering towards the observer. As with the EBT, 
we must also account for the probability of any scat¬ 
tered/emitted radiation reaching the observer via a suitable 
mean exponential factor, < >. Thus, in the TBT, we 

compute synthetic observables via 



w_^_ 

At Av 



< e 



(19) 


where the summation is again over i’—packets in the selected 
frequency and time bins. Formally, the exponential weight 
factor of a r’—packet in the TBT should be computed as 

i-\-/\i 

f ^ dl' 

< >= - -- , ( 20 ) 

where the integral runs over the trajectory path Al. To hrst 
order, however, this can be approximated by generating the 
i’—packet at the mid point of Al and computing the expo¬ 
nential factor from its hight (as described in Section 2 . 2 . 2 ); 

i.e. 


^ g — '^esc g —Tesc(^ + ^V2) (21) 

In principle, dP/dQ\TBT can be formulated to account for all 
(effective) scattering/huorescence processes. However, since 
we are primarily interested in studying contributions to the 
emergent polarisation spectrum, we focus only on electron 
scattering, for which 


d^ 

dTt 


= j- P{e,ii) hcAl 

47r 


( 22 ) 


where ksc is the scattering attenuation coefficient. 

A key difference between the EBT and the TBT is that, 
in the latter, every trajectory element of the r—packets con¬ 
tributes to the synthetic observables, whereas, in the former, 
only physical interaction events contribute. For instance, 
in the limit of optically thin ejecta many more i’—packets 
would contribute to the emergent spectrum in the TBT com¬ 
pared to the EBT (see Fig. 3). However, a drawback of the 
TBT is that resc{l + A//2) should describe the mean proba¬ 
bility of escape for points along the trajectory element (see 
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Figure 4. The geometry adopted by Hillier (1994) and used for 
our test calculations. Monte Carlo quanta are created inside a 
spherical shell of radius i^min, propagate into a region with a pro¬ 
late density distribution and are free to escape when they reach 
the outer spherical shell at i^max = SO.ORmin- -^e(r,/5) is the 
electron density distribution and Ago = Ae(Amin, 0). 


equation 21 ), rather than the exact probability of escape 
from the interaction point, as in the EBT. For an r—packet 
trajectory with moderate optical depth (r > 1 ), this ap¬ 
proximation may lead to a poor estimate of the observables. 
Breaking the trajectory into smaller paths (with At < 1) 
and generating packets at each midpoint may be required, 
slowing down the code and making the TBT less efficient 
then the EBT (see Section 3.2). 


3 TEST CODE 

In the following we present a simple test code, with the aim 
of validating our polarisation scheme (Section 3.1) and se¬ 
lecting the most suitable of the techniques described in Sec¬ 
tion 2.2 to synthesise spectra with low Monte Carlo noise 
(Section 3.2). In this code, packets are generated with null 
polarisation in a small inner sphere (to mimic a point source) 
and then allowed to propagate into an envelope where either 
interactions with electrons or continuum absorptions can oc¬ 
cur. Here time evolution for the ejecta is neglected. 


3.1 Polarisation scheme validation 

To validate the polarisation scheme, we hrst focus our at¬ 
tention on the DCT and choose to reproduce a simple con- 
hguration described by Hillier (1994). As shown in Fig. 4, a 
point source is surrounded by a detached spherical shell with 
inner radius Amin = 2.0 and outer radius Amax = 30.0Amin, 
with a prolate density distribution Ne{r, (3) such that 

aeNe{r,/3) = (1 + 10 cos^^) , (23) 

where cTg is the Thomson cross section and r and /3 express 
the radius and the polar angle inside the envelope. The xo 
parameter is related to the solid-angle averaged (from in¬ 
ner to outer boundary) optical depth, Tave = 2 . 888 x 0 , and 
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Figure 5. Continuum polarisation as a function of xo and r^ve for 
4 different viewing angles for the setup described in Section 3.1. 
Symbols indicate predictions from our test calculations, while the 
lines are reported from Hillier (1994) for comparison. The Monte 
Carlo noise error bars are not shown since they are smaller than 
the symbol sizes. 


Figure 6. Continuum polarisation as a function of the albedo 
(black points) for the setup described in Section 3.1, together 
with the predicted curve (dashed line) from Hillier (1994). Here 
the observer’s inclination is 90° and xo = 0. The Monte Carlo 
noise error bars are not shown since they are smaller than the 
symbol sizes. 


can be varied to investigate the impact of different scatter¬ 
ing optical depth on the polarisation signal. Neglecting ab¬ 
sorption and assuming a pure electron scattering envelope, 
the continuum polarisation as a function of xo is shown in 
Fig. 5 for four different viewing angles i (22.5°, 45°, 67.5° 
and 90°). The agreement between our predicted values and 
the expected curves from Hillier (1994) is encouraging. We 
also carried out calculations in which we include continuum 
absorption opacity. These show good agreement with the 
predicted dependence of the continuum polarisation on the 
albedo (ratio of the scattering to the total opacity, see Fig. 
6 ). 

3.2 Comparison between different techniques 

A convenient means to compare the three techniques for ex¬ 
tracting observables outlined in Section 2.2 is by studying 
their accuracy in reproducing continuum polarisation for a 
given configuration. We did this by repeatedly running our 
test code a number (iVsim = 500) of times for each tech¬ 
nique (with different random number seeds determined by 
the wall-clock time) and comparing the distributions of po¬ 
larisation values obtained using each method. For these ex¬ 
periments, we chose a configuration in which a point source 
illuminates a surrounding atmosphere, chosen to be a con¬ 
stant density oblate ellipsoid with axis ratio of two, i.e. 

^ + S + 4 = l, a = b = 2c . (24) 

In all the simulations, 10® packets have been created, a pure 
scattering atmosphere has been assumed (i.e. no line opacity 
and albedo equal to 1) and the scattering coefficient has been 
set to ksc = 1/a and kept constant throughout the ejecta. 
As mentioned in Section 2.2.1, continuum polarisation levels 
estimated from the DCT may be either inaccurate or have 


large uncertainties depending on whether the number of an¬ 
gular bins is small or large, respectively. We found that a 
value of 51 angular bins provides a reasonable compromise 
and we therefore adopt this in the DCT calculations. The 
results of the comparison between the three techniques are 
shown in Fig. 7 and reported in Table 1. 

For an observer along the y-axis (Fig. 7, left panel), the 
projection of the model along the line-of-sight is circular and 
we expect to find null polarisation (see above). The Q and 
U values in every technique are indeed consistent with zero 
but, because of Monte Carlo noise, a distribution of values 
is obtained. The width of this distribution provides a conve¬ 
nient quantification of the Monte Carlo noise introduced by 
each method. As expected, the standard deviations obtained 
with the EBT and TBT methods are smaller than with the 
DCT method, by factors of ~ 6.6 and ~ 8.6, respectively. 
Given that Monte Carlo noise is expected to scale with the 
square root of the number of packets, the DCT could reach 
the same signal-to-noise ratio as the EBT (TBT) with a 
factor of ~ 45 (~ 75) more packets. Although the i’—packet 
routine introduces computational overheads in the EBT and 
TBT (see Table 1), the direct counting approach would still 
be a factor of ~ 7.5 (~ 2.5) slower than the i’—packets tech¬ 
nique. We also note that, even with the same signal-to-noise 
ratio, the DCT would be less accurate in predicting polar¬ 
isation values because of the need to average contributions 
from different viewing angles (via angular binning): indeed, 
closer inspection shows that the Q distribution for the DCT 
is shifted towards positive values (in the binning approach 
contributions to the spectra come from r—packets escaping 
close to, but not exactly along, the z direction and thus the 
average value of Q is slightly positive). 

If viewed down the x-axis (right panel), the projection 
becomes an ellipse and thus we expect to see a polarisation 
signal. The three techniques agree in reproducing a contin¬ 
uum polarisation of p ~ 3.8 per cent but, again, the DCT 
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Figure 7. Q and U continuum polarisation distributions of Nsim = 500 runs with the DCT (top panels), the EBT (middle panels) and 
the TBT (lower panels). The adopted geometry is an oblate ellipsoid with axis ratio of two, as described in Section 3.2. The system is 
viewed along the minor axis (circular projection, left panel) and the major axis (elliptical projection, right panel). For each distribution, 
a solid vertical line indicates the average value, and the grey shaded area marks ± one standard deviation, cr. The average values and 
the standard deviations of each plot are reported in Table 1 together with the average runtimes, t. 


Table 1. Average values and standard deviations of the distribu¬ 
tions of Q and U values predicted on 500 simulations by the DCT, 
EBT and TBT. The system is an oblate ellipsoid with axis ratio 
of two viewed along the 2 ;-axis (circular projection) and along the 
x-axis (elliptical projection). The averaged runtime, t, is reported 
for each distribution. 


Circular 

Q ± fjQ (per cent) 

U ± ajj (per cent) 

t (s) 

DCT 

0.16 ± 0.23 

0.00 ± 0.24 

3.9 

EBT 

0.00 ± 0.04 

0.00 ± 0.04 

22.3 

TBT 

0.00 ± 0.03 

0.00 ± 0.03 

123.0 

Elliptical 

Q ± fjQ (per cent) 

U ± ajj (per cent) 

t (s) 

DCT 

3.81 ± 0.28 

-0.02 ± 0.29 

3.9 

EBT 

3.82 ± 0.04 

0.00 ± 0.04 

26.3 

TBT 

3.81 ± 0.03 

0.00 ± 0.03 

144.5 


gives a much broader distribution of values, indicating that 
it is more severely affected by noise in the simulation. In or¬ 
der to reach the same standard error of the mean of the EBT 
(TBT), the DCT would require a factor of ~ 50 (~ 100) 
more packets, making the total runtime a factor of ~ 7.5 
(~ 3) longer than the r'—packets scheme. 

This simple comparison shows that the packet ap¬ 
proaches are more precise in estimating polarisation, allow¬ 
ing us to reach a given Monte Carlo noise with many fewer 
Monte Carlo quanta (and substantially shorter runtimes) 
than the simple DCT would require. As already anticipated 
in Section 2.2.3, the EBT is indeed more efficient than the 
TBT because the runtime for the latter is limited by the 
need of breaking r—packet trajectories with moderate opti¬ 
cal depth (r > 1) into smaller paths (with At <C 1), in order 
to give accurate results for Tesc- We note that, although we 
have carried out polarisation tests here, this improvement 
in Monte Carlo noise could also be exploited for extracting 
high-quality observables of any sort. 


4 ARTIS 

In this Section we describe the implementation of our polar¬ 
isation scheme into the three-dimensional, time-dependent 
radiative transfer code ARTIS (Kromer & Sim 2009) and 
test it for one-dimensional and two-dimensional models. Sec¬ 
tion 4.1 outlines the implementation of the different tech¬ 
niques described in Section 2.2. In Section 4.2 we test 
the code using the one-dimensional W7 explosion model 
(Nomoto et al. 1984; Iwamoto et al. 1999) and check the ac¬ 
curacy of the t’—packet technique in computing spectra and 
reproducing continuum polarisation consistent with zero. Fi¬ 
nally, in Section 4.3 we apply the new version of the code 
to two-dimensional ellipsoidal models to investigate the im¬ 
pact of simple aspherical geometries on line and continuum 
polarisation for different viewing angles, and compare our 
results to those of similar studies made using other codes. 

4.1 Implementation 

Polarisation is implemented into ARTIS by assigning a 
Stokes vector to each r—packet and by transforming this 
according to the physical process the packet undergoes (see 
Section 2.1). For the DCT, the same binning approach al¬ 
ready used in the code for spectra and light curves is ex¬ 
tended to compute polarisation spectra. 

As mentioned in Section 2.2.3, the r'—packet TBT 
should yield spectra with lower Monte Carlo noise compared 
to the EBT since contributions to the polarisation spectra 
come from every event in the r—packet histories, including 
both physical interactions and numerical events (e.g. pack¬ 
ets crossing grid cell boundaries). However, accurate results 
from the TBT require that care is taken in the calculation of 
Tesc, which can introduce a large computational overhead for 
complicated opacity distributions. Indeed, our test calcula¬ 
tions (Section 3.2), suggested that this additional computa¬ 
tional overhead can ultimately make the TBT less efficient 
than the EBT. Consequently, here we have chosen to imple- 
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ment a i’—packet routine using the EBT that can be used 
to compute synthetic observables from ARTIS. 

The u—packet routine allows us to compute flux and 
polarisation spectra for multiple observers, simply by using 
a loop to generate u—packets over a set of different viewing 
angles. Several input parameters can be chosen to optimise 
performance in the calculations. First, the calculation of the 
optical depth Tesc, see equation (17), can be stopped when 
the u—packet reaches a maximum value u— packets 

with high optical depth to the boundary would make van¬ 
ishingly small contributions to the final spectrum (because 
of the exponential factor, see equation 15) and can thus 
be neglected. Since the computation cost of the u—packet 
methods is dominated by the calculation of Tesc, the run¬ 
time is strongly affected by this parameter. We have carried 
out test calculations to verify that setting = 10 and 

neglecting u—packets with higher optical depth does not af¬ 
fect the hnal result, and adopted this as our default value 
for all the calculations presented here. Our implementation 
also includes the option to generate u—packets only in a 
selected spectral interval. Because spectropolarimetric ob¬ 
servations usually cover the optical region of the spectrum, 
our default wavelength range is 3500 — 10 000 A. Given that 
much of the runtime of the code can be consumed in comput¬ 
ing Tesc for packets in the bluer regions, which easily reach 
because of strong iron-line blanketing, this particular 
cut in wavelength speeds up the calculations by a factor of 
~ 4 compared to calculations with no wavelength cut. Fi¬ 
nally, the u—packet routine can be switched on or off for 
timesteps as chosen by the user (note that the activation 
or deactivation of u—packets has no effect on the r—packet 
propagation). 

4.2 W7 model 

Although much effort has been recently directed at devel¬ 
oping multi-dimensional explosion models (Rosswog et al. 
2009; Jordan IV et al. 2012; Kushnir et al. 2013; Moll et al. 
2014; Fink et al. 2014), the one-dimensional parameterised 
deflagration model W7 is still widely used since its compo¬ 
sition and structure provide reasonable agreement with ob¬ 
servations of (“normal”) SNe la (Kasen et al. 2006; Kromer 
V Sim 2009; Jack et al. 2011; van Rossum 2012; Gall et al. 
2012). Here we calculate ARTIS flux and polarisation spec¬ 
tra for this model aiming to: (i) compare the accuracy of the 
DGT and the EBT in reproducing flux spectra at different 
epochs; (ii) test our polarisation implementation on a spher¬ 
ically symmetric system for which null polarisation spectra 
are expected. 

For this calculation we simulate 8 x 10^ Monte Garlo 
packets and compute spectra over 111 logarithmic time- 
steps from 2 to 120 days after explosion. We bin the final 
spectra of both the DGT and the EBT in logarithmic wave¬ 
length bins with ^ = 3.912 x 10“^. For this test, local 
thermodynamic equilibrium (LTE) has been assumed for all 
time-steps^. The calculation was carried out by mapping the 

^ We note that LTE is a crude approximation, especially for 
epochs after maximum light. We will, however, confine most of 
our discussion to relatively early epochs when the LTE approxi¬ 
mation should be reasonable. 



Figure 8. Spectra for the W7 model at 15, 20 and 25 days 
after explosion computed with the DOT (black lines) and 
EBT (red lines). The spectra are computed for an observer at 
"^obs = (0, 0, !)• The model supernova is assumed to be at 1 Mpc. 

spherically symmetric W7 model onto a 100^ Gartesian grid, 
through which the packets were propagated. The u—packet 
EBT is activated from 10 to 30 days after explosion and only 
for r—packets with emergent rf wavelength between 3500 
and 10 000 A. Spectra for the EBT are computed for the 
viewing angle riobs = (0,0,1), although we note that (since 
the model is spherically symmetric) the choice of observer 
orientation here is arbitrary. Gompared to the DGT, the 
runtime penalty associated with using the u—packet routine 
in the EBT is found to be less than a factor of two, with the 
advantage that the number of packets contributing to the 
emergent spectrum is a factor of ~ 115 higher. 

Fig. 8 shows spectra calculated with the EBT at 15, 20 
(around B band maximum light) and 25 days after explo¬ 
sion. Angle-resolved (10 angle bins^) spectra obtained with 
the DGT are shown for comparison. We note that the cal¬ 
culation of the direct counting flux spectrum is exactly the 
same as in the previous version of ARTIS (Kromer & Sim 
2009), with the exception that electron scattering is now 
treated fully via the scattering matrix in equation (9) rather 
than assuming isotropic scattering. The agreement between 
the two techniques is very good, with the EBT being much 
less affected by Monte Garlo noise, as expected. To estimate 

In Section 3.2 we chose a number of 51 viewing-angle bins as 
a reasonable value to obtain accurate angle-resolved results with 
relatively low Monte Carlo noise for a simple ellipsoidal configura¬ 
tion. However, the number of packets used here requires a smaller 
number of bins in order to achieve a reasonable level of Monte 
Carlo noise in the spectra. Reducing the number of viewing-angle 
bins to 10 does not affect the accuracy of the results (given that 
the observables in a ID model are the same for different viewing 
angles) but instead merely decreases the Monte Carlo noise in the 
predicted spectra. 
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Figure 9. Accuracy of the DCT and the EBT in reproducing 
continuum polarisation consistent with zero. The W7 model has 
been used for this test calculation. Polarisation spectra are com¬ 
puted at 20 days after explosion with the DCT (black lines) and 
EBT (red lines). The increase in Monte Carlo noise at longer 
wavelengths is due to the lower flux in the spectrum (see Fig. 8). 


the Monte Carlo noise in the spectra, we use the fact that 
the calculation has been carried out on multiple cores which 
provides us with a set of independent estimates for any given 
observable. In particular, we divide the simulation outputs 
into eight subsets, each comprising one eighth of the cores, 
and calculate an emerging spectrum for each of them. Spec¬ 
tral differences between different subsets are representative 
of the Monte Carlo noise and estimated by computing resid¬ 
uals from a mean spectrum. The standard deviation of the 
residual in the packet spectrum is 13.3 times smaller than 
that calculated for the angle-resolved direct counting spec¬ 
trum. 

Polarisation spectra around maximum light in the B 
band are reported in Fig. 9. As expected from a one¬ 
dimensional model, the average Q and U throughout the 
whole spectral range are consistent with zero for both tech¬ 
niques, with the signal-to-noise decreasing towards the red 
because of the lower flux level. The decrease in Monte Carlo 
noise when comparing the DCT to the EBT is remarkable: 
the standard deviation in the Q (U) spectrum is a factor of 
14.1 (13.7) larger in the former compared to the latter, in 
good agreement with our hndings from the flux spectra. This 
simple comparison clearly shows that the r'—packet tech¬ 
nique is superior for producing accurate polarisation levels. 

We note that the factor by which the noise improves 
does depend on the number of angular bins used (the im¬ 
provement is less dramatic - but still signihcant - if fewer 
bins are used). However, our choice of 10 bins is rather con¬ 
servative (c.f. Section 3.2 where it was found that ~ 51 bins 
were required for accurate representation of a simple 2D 
model). 


Figure 10. Spectra for the PEM (top panel) and the OEM 
(bottom panel) calculated with the EBT at 19 days after explo¬ 
sion. Black/green lines are for an observer orientation along zjx 
iP'ohs.il '^obs ,2 in the text). Scaled projected surfaces are shown 
for each viewing angle. 


4.3 Ellipsoidal toy model 


In this section, we follow previous studies (e.g. Hoflich 1991; 
Kromer & Sim 2009; Dessart & Hillier 2011) and use el¬ 
lipsoidal models as a starting point to explore aspherical 
geometries. We use a model equivalent to that of Kromer & 
Sim (2009), which has a prolate ejecta morphology (PEM) 
and also consider a similar model with an oblate ejecta mor¬ 
phology (OEM). Specihcahy, we assume ellipsoidal isoden¬ 
sity surfaces with density prohle 


p(0 « 


®^p(-^) 

0 


max 

max 


(25) 


The parameter ^ is dehned in terms of the components of 
velocity in cylindrical polar coordinates v — as 



OEM 

PEM 


(26) 


where = 2750 km s“^ and h is the ratio between the semi¬ 
major and semi-minor axis. Here we limit our study to an 
axis ratio h — 2 and hx the maximum velocity parameter 
^max = 13 750 km s“^. We adopt a composition for both 
models that is roughly appropriate for SNe la. Specihcahy, 
the total mass and chemical yields of the ejecta are chosen to 
be the same as for the W7 model and a stratihed composition 
with three ellipsoidal zones {h = 2) is assumed. The model is 
set up by hhing the ejecta from the centre to accommodate 
the W7 yields of the different element groups: the innermost 
region is hhed with “iron group” elements (21 ^ Z ^ 30), 
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Figure 11. Flux and polarisation spectra at 19 days after explosion for the PEM (left panels) and the OEM (right panels). The observer 
is placed at nobs,2 (along x). Identification between polarisation features and lines in the spectrum are shown with blue vertical lines. 
Scaled projected surfaces are shown in green. 
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Figure 12. Colour maps of normalised / (left panels), Q (middle panels) and U (right panels) distributions projected on the velocity 
plane (vy^Vz). The maps are computed for the PEM using the EBT and selecting the emergent n—packets between 16.5 and 21.5 days 
after explosion and in the wavelength regions 3500 — 6000 A (upper panels) and 6400 — 7200 A (lower panels). Solid lines mark the outer 
boundary of the iron-group-element zone (inner white ellipse) and the maximum velocity parameter ^max = 13 750 km s“^ (outer black 
ellipse). The intensity distribution in the blue is more circular than the projected density contour: Q is dominated by contributions along 
the minor axis, leading to a positive polarisation level. In contrast, the intensity distribution in the red is more similar to the projected 
density contour: Q is dominated by contributions along the major axis and therefore biased towards a negative value. 
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Figure 13. Evolution of the relative contribution of different 
scattering processes to the observed spectrum. Fractions are cal¬ 
culated for the PEM (left panels) and the OEM (right panels) 
with the DOT by selecting escaping packets based on their last 
interaction(s) prior to escape. The fraction of packets that under¬ 
went a depolarising interaction process (bound-bound, bound-free 
or free-free emission) as last interaction is shown in red. The con¬ 
tribution from packets that had a single electron scattering inter¬ 
action since their last depolarising interaction is indicated in blue, 
and packets that suffered multiple electron scattering events prior 
to escaping are show in black. Upper panels show contribution in 
the spectral region 3500 — 6000 A, lower panels in the wavelength 
range between 6400 and 7200 A. 

the middle with intermediate-mass elements (9 ^ ^ 20) 

and the outermost with low-mass elements (Z ^ 8). The 
transitions between the different layers are at ^ ~ 8000 and 
^ ~ 10 500 km s“^. Relative abundances of the elements 
inside each zone are kept fixed to the W7 values. 

LTE radiative transfer calculations have been per¬ 
formed over 111 logarithmic time-steps from 2 to 120 days 
after explosion. 2.4 x 10^ and 4 x 10^ Monte Carlo quanta 
have been generated for the PEM and OEM, respectively. 
Since the redder regions of the polarisation spectra are typ¬ 
ically noisier due to the lower flux (i.e. fewer Monte Carlo 
quanta per frequency bin; see Fig 9), an additional simula¬ 
tion has been carried out for the PEM (OEM), with 8 x 10^ 
(4 X 10^) Monte Carlo quanta and with the EBT routine 
called only for A > 6000 A^. The two simulations have 
thus been combined to produce final spectra for the EBT 
in the whole range between 3500 and 10 000 A. Spectra are 
computed with the EBT from 10 to 30 days after explo¬ 
sion and for two extreme viewing angles: along the z-axis, 
^obs,i = (0,0,1), and along the x-axis, nobs ,2 = (1,0,0). 

4-3.1 Flux and polarisation spectra 

In Fig. 10 we compare the i’—packet total flux spectra at 
19 days after explosion for the two ellipsoidal models. We 
find the same strong viewing-angle dependencies reported 

^ This cut in wavelength considerably speeds up the calculation 
since the n—packet routine is called a factor of ~ 25 times fewer 
compared to calculations with the entire range 3500 — 10 000 A. 


Figure 14. Flux spectrum (solid black line) and Q polarisation 
spectrum (red line) around the Si ll A5979 and Si ll A6355 features 
for the PEM viewed along the a:-axis. Rest wavelengths of the two 
lines are marked by vertical dashed lines. Inverted P-Cygni pro¬ 
files for the two silicon lines can be identified in the Q spectrum. 

by Kromer & Sim (2009). For a given morphology, packets 
escaping along the major axis see a velocity twice as large 
compared to the minor axis and the corresponding spectrum 
is therefore characterised by broader features and stronger 
line blending; moreover, the spectrum viewed along the ma¬ 
jor axis is fainter since the projected area along this axis 
is smaller and the typical opacity is higher. The same ge¬ 
ometrical arguments can also be used to compare spectra 
for the two different geometries: spectra viewed down the 
minor (major) axis are qualitatively similar, because pack¬ 
ets see the same velocity range, but the prolate ellipsoid is 
fainter than the oblate due to the smaller projected surface. 

Polarisation spectra for the observer orientation riobs,! 
are consistent with zero for both models, reflecting the over¬ 
all spherical symmetry of the projected surface. As shown in 
Fig. 11, however, observer orientations from which the model 
has an elliptical projected surface produce a clear polarisa¬ 
tion signal in Q. U remains consistent with zero because the 
model is axi-symmetric, and the calculated U spectrum can 
be used as a convenient proxy for the Monte Carlo noise in 
the Q spectrum. 

Sign reversals from shorter to longer wavelengths are 
found in the Q spectrum for both the PEM and the OEM, a 
behaviour that can not be explained by the simple picture of 
an optically thin electron scattering atmosphere illuminated 
by a point source. In the latter, one would expect the overall 
polarisation to be negative (positive) for the PEM (OEM), 
with a polarisation decrease across the lines because of flux 
dilution. Instead, as found by previous studies (Dessart & 
Hillier 2011; Pat at et al. 2012), the results of full calculations 
are more complex and sign reversal in polarisation spectra 
can arise. These complexities can be ascribed to variations 
in thermalisation depth with wavelength (see below for ex¬ 
planation) and highlight the need for realistic calculations 
beyond simple toy atmosphere geometries for the interpre¬ 
tation of data. 

Fig. 12 shows the intensity and polarisation distribu¬ 
tions projected on the velocity plane {vy,Vz). The maps 
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Figure 15. Flux and polarisation spectra for the PEM (left panels) and the OEM (right panels) calculated for a viewing angle nobs,2 
(along x) at 14 (orange), 19 (red) and 24 (blue) days after explosion. Scaled projected surfaces are shown in green. 


have been calculated for the PEM selecting the emergent 
r’—packets between 16.5 and 21.5 days after explosion and 
in the spectral regions 3500 — 6000 A and 6400 — 7200 A. 
In both wavelength intervals, the intensity emission region 
in projection is less elliptical than the density contour, and 
this is a stronger effect in the blue. This behaviour can be 
ascribed to the relative contributions of the line and the elec¬ 
tron scattering opacities in different regions of the spectrum 
(see Fig. 13): the blue region is dominated by line opaci¬ 
ties and thus the intensity distribution in projection is more 
circular than elliptical®; in contrast, the red region is free 
from strong line opacities (see also Pinto & Eastman 2000, 
Kasen et al. 2004 and Patat et al. 2009) and therefore the 
projected intensity distribution is more similar to the ellip¬ 
tical density contour. Because contributions to Q are typi¬ 
cally positive from regions along the minor axis and negative 
along the major axis, different distributions in intensity lead 
to different values of the overall Q polarisation: in the blue, 
polarisation in Q is dominated by contributions along the 
minor axis and thus biased towards a positive value, while 
in the red the Q emission is stronger along the major axis 
and thus biased toward a negative value. The same argu¬ 
ments explain why the Q polarisation in the OEM changes 
from negative values in the blue to positive in the red. 

The blue region is also characterised by strong 
(~ 1 — 2 per cent) polarisation across spectral lines: polari¬ 
sation peaks are associated with the blue-shifted absorption 


® This is because, among all the packets created at a given iso¬ 
density surface, those at highest projected velocities (i.e. around 
the major axis of the ellipsoid) sweep out the largest velocity 
range on their journey to the observer, and therefore encounter 
the greatest line opacity. 


trough, where contributions from the weakly polarised cen¬ 
tral source are scattered out of the line of sight by the line. 
In contrast, a decrease in polarisation is found in the emis¬ 
sion wing of the P-Cygni profile, where line scattering brings 
extra unpolarised packets into the line of sight. This leads 
to the inverted P-Cygni profile shape in the Q spectrum, as 
discussed by Jeffery (1989). This is clearly visible in the two 
Sill lines at ~ 5979 and ~ 6355 A (see Fig. 14). 

4-3.2 Spectral evolution and light curves 

The spectra of both ellipsoidal models are shown in Fig. 15 
for three epochs (14, 19 and 24 days after explosion). The in¬ 
tegrated luminosity in the wavelength range 3500 — 10 000 A 
peaks at ~ 19 days after explosion in the PEM, whereas the 
maximum is reached earlier in the OEM (~ 14 days after 
explosion, see Fig. 16). To quantify the time evolution of 
the polarisation, we have also computed polarisation light 
curves Q*(t) and f/*(t) by integrating Q and U values over 
chosen wavelength regions (Ai to A 2 ): 

Q*(t) = f ' Q(A,t) dA , f/*(t) = /* ' U{X,t) dX . (27) 

Jxi Jxi 

As expected, U* remains consistent with zero at all times. If 
we consider the entire synthetic spectrum (i.e. Ai = 3500 A, 
A 2 = 10 000 A), Q* in the PEM (OEM) evolves from nega¬ 
tive to positive (positive to negative) values as we go from 
early to late epochs (see Fig. 16). This behaviour is due to 
changes in the relative contributions of the blue and red re¬ 
gion as a function of time, and can easily be understood 
from polarisation light curves computed for the spectral in¬ 
tervals between Ai = 3500 A and A 2 = 6000 A and between 
Al = 6400 A and A 2 = 7200 A (see Fig. 16). As can be 
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PEM OEM 






Figure 16. Polarisation light curves between 10 and 30 days 
after explosion for the PEM (left panels) and the OEM (right 
panels) viewed along the ai-axis. Only Q* is reported here since 
U* is consistent with zero for both models. Grey lines represent 
contributions from the whole spectral range for which i;—packets 
were calculated (3500 — 10 000 A), whereas blue and red lines 
are for packets escaping at short (3500 — 6000 A) or longer 
(6400 — 7200 A) wavelengths. Spectral flux integrated in the whole 
range is reported in the lower panels. One-sigma Monte Carlo 
noise error bars are derived using the procedure outlined in Sec¬ 
tion 4.2. Some error bars are not visible because they are smaller 
than the symbol sizes. 


Wavelength (A) 

Figure 17. Q polarisation spectrum for the PEM (in red) calcu¬ 
lated ~ 4 days after B band maximum light. For comparison the 
black line shows the polarisation spectrum of SN 2004dt (Leonard 
et al. 2005) at the same epoch. Given that the PEM is axisym- 
metric, the polarisation spectrum of SN 2004dt calculated along 
the dominant axes, is shown here. Polarisation levels predicted 
by the PEM across the lines are comparable (within a factor of 
two) to those observed, while the polarisation in the continuum 
(grey shaded area) is too high. 


anticipated from the polarisation spectra, Q* is positive in 
the blue and negative in the red for the PEM, whereas the 
opposite is true in the OEM. The OEM light curve evolves 
more rapidly than the PEM, having reached peak flux at 
around 14 days and then starting to decline. This more 
rapid evolution is also clearly evident in the degree of po¬ 
larisation, which changes much more noticeably over this 
time period for the OEM. In particular, polarisation in the 
pseudo-continuum region between 6400 and 7200 A is ap¬ 
proximately constant (Q*//* ~ — 1 per cent) in the PEM, 
whereas significant evolution is found for the OEM. 

Although the ellipsoidal models studied here are not 
very realistic, they do qualitatively reproduce the main fea¬ 
tures of SN la polarisation spectra. As shown in Eig. 17, the 
PEM predicts an overall small polarisation signal through¬ 
out the spectrum {p <1.5 per cent) and polarisation levels 
across the lines comparable (within a factor of two) to those 
observed in SN la. Of course, the comparison with data is far 
from perfect: the polarisation predicted in the continuum is 
too high (because of the strong asymmetry in the toy model) 
and the velocities of the line features are too small. Such dis¬ 
crepancies come as no surprise, given the simplicity of the 
model. In future studies, we will make quantitative com¬ 
parisons to data with results from polarisation calculations 
performed for real explosion models. 


5 CONCLUSIONS 

We have described a technique for modelling polarisation for 
multi-dimensional supernova explosion simulations, and im¬ 


plemented it in the radiative transfer code ARTIS (Kromer 
& Sim 2009). Our method uses an approach inspired by 
Lucy (2005), and related to those used by Long & Knigge 
(2002), Sim et al. (2010) and Kerzendorf & Sim (2014), for 
extracting observables: viewing-angle spectra are obtained 
by summing contributions from i’—packets generated at each 
r—packet interaction point and escaping to a chosen ob¬ 
server orientation (Event-Based Technique, EBT). These es¬ 
caping i’—packets can be used to construct synthetic observ¬ 
ables (spectra, light curves, polarisation spectra) that have 
substantially reduced Monte Carlo noise compared to spec¬ 
tra obtained by direct binning of escaping r—packets. We 
also investigated a higher-order Monte Carlo noise reduction 
approach, based not on r—packet interaction sites but on 
r—packet trajectory elements (Trajectory-Based Technique, 
TBT). Initial results, however, suggest that this approach 
may be less computationally efhcient. 

We validated our polarisation scheme using an idealised 
test code in a simple configuration, and found good agree¬ 
ment with predictions from Hillier (1994). Applying the 
same idealised test code to a simple ellipsoidal toy model, we 
then verified that continuum polarisation levels calculated 
with the EBT agree with values predicted by direct binning 
of the escaping r—packets. We implemented the EBT in AR¬ 
TIS and tested it for a model with a realistic SN la compo¬ 
sition and opacity (the spherically symmetric W7 model): 
as expected, the r'—packet method could accurately repro¬ 
duce the synthetic spectrum obtained by direct binning of 
emergent Monte Carlo quanta and also predict polarisation 
consistent with zero. However, the EBT is much less affected 
by Monte Carlo noise (with typical signal-to-noise ratios a 
factor of ~ 13 higher than those obtained with the direct 
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binning approach) and thus more suitable to reproduce very 
weak signals (e.g. polarisation levels observed in SNe la). 

Finally, we synthesised flux and polarisation optical 
spectra with the EBT for prolate and oblate ellipsoids with 
axis ratio of two, using typical SNe la velocities and com¬ 
positions (including composition layering). As expected, we 
obtained null polarisation spectra when the projected sur¬ 
face on the plane of the sky is circular. In contrast, as- 
pherical projected areas yield a polarisation signal (typically 
~ 1 per cent) in both morphologies. The polarisation is char¬ 
acterised by sign reversals across the spectrum and peaks 
associated with troughs of strong optical features. This be¬ 
haviour is consistent with results of previous studies using 
similar ejecta morphologies (Hoflich 1991; Dessart & Hillier 
2011; Patat et al. 2012) and is ascribed to variations in ther- 
malisation depth with wavelength. At the epochs we studied 
(14 — 24 days post-explosion), the evolution of polarisation 
spectra is more dramatic for the oblate than the prolate mor¬ 
phology, both in the continuum and in the line polarisation 
levels. 

In this paper we have focused on developing our tech¬ 
nique and testing its accuracy in calculating intensity and 
polarisation spectra for one- and two-dimensional models. 
This study has laid the groundwork for future calculations 
in which we will exploit the multi-dimensional capability of 
ARTIS and calculate polarisation spectra for a set of con¬ 
temporary SN la explosion models. Such calculations will 
help to identify geometric discriminants between models and 
to make comparisons between their predictions and spec- 
tropolarimetric data more reliable. 
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